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ABSTRACT 


The fitting of radial velocity curves is a frequent procedure in binary stars 
and exoplanet research. In the majority of cases the fitting routines need to be 
fed with a set of initial parameter values and priors from which to begin the 
computations and their results can be affected by local minima. We present 
a new code, the rvfit code, for htting radial velocities of stellar binaries and 
exoplanets using an Adaptive Simulated Annealing (ASA) global minimization 
method, which fastly converges to a global solution minimum without the need 
to provide preliminary parameter values. We show the performance of the code 
using both synthetic and real data sets: double-lined binaries, single-lined bi¬ 
naries, and exoplanet systems. In all examples the keplerian orbital parameters 
htted by the rvfit code and their computed uncertainties are compared with 
literature solutions. Finally, we provide the source code with a working example 
and a detailed description on how to use it. 

Subject headings: Data Analysis and Techniques — Stars — Extrasolar Planets 
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Introduction 


Precise, absolute masses of stars and planets can be only measured if they are part 
of binary (star-star or star-planet) systems. The masses are derived by fitting the radial 
velocity (RV) curves of those systems to the keplerian orbital equations. Those equations 
need to be solved numerically using multi-parameter minimization techniques. 

The fitting of radial velocity curves is an a-priori straightforward procedure that gets 
complicated by the need to explore a wide multi-parameter space and by the existence 
of many potential local minima, which can yield to incorrect solutions. Local minima 
are a well known issue in mathematical optimization problems and a lot of work has 
been done over the years to overcome this limitation when trying to £t multi-parameter 


functions. Simultaneous, mul ti-parameter minimization techniques, such as 


Marquardt (LM) algorithm flLevenberg 


1944J: 


Marquardt 


1963 


Press et ah 


he L evenberg- 


19921) . also 


known as t 


le damped leas 


^-squares (DLS) m ethod, or the Nelder & Mead Simplex (NMS) 


algorithm flNelder &: Mead 


1965 


Murtv 


19831) . are widely used in multi-parameter function 


optimization problems in many science fields. 

Different minimization methods are implemented in the most widely used binary and 
exoplanet modelling softw ares. Some of those softwares are the Wilson-Devinney code 


(jWilson fc Devinnev 


197ll) which is the de facto standard code to analyze eclipsing binary 


stars and implements the Differential Corrections (DC) method which is derivative based 
and is used to look for local conver gence. A wrapper for th e Wilson-Devinney Method, 


WD2007, was built by J. Kallrath ( Milone fc Kallrath 


Anne aling (SA) in the Boltzmann version (see e. g 


20081) and implements Simulated 


Metropolis et a 


19831 ). PHOEBE (PHysics Of Eclipsing BinariEs; 


Prsa fc Zwitter 


1953 


Kirkpatrick et ah 


20051) is a code which 


models light and RV curves. PHOEBE relies in the Wilson-Devinney code and acts as a 
front-end interface adding new capabilities. Another commonly used software is Nightfall 
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flWinhmannl l201ll) . which is a code to build synthetic light and RV curves of eclipsing 


binary stars and can fit a model to photometric and RV data. Ni ghtfall uses 


Algorithm’ for local optimizacion and SA for global optimization flWichmann 


he ’S implex 


1999|). 


In recent years, the discovery of hundreds of exoplanets using RV techniques has lead 
to the development of several software packages to fit RV and transit photometric curve s. 


One example is the Systemic Console flMeschiari et al 


2009 


Meschiari fc Tjanghlin 


201011 . 


which focuses in the htting of exoplanet RV and transit curves and includes the LM, NMS 
and SA algorithms. 

All those software packages include methods for local and global minimisation, SA 
among them. SA is a multi-parameter mi nimization technique th at draws on metallurgic 


optimal cooling processing methods (see 


Kirkpatrick et al. 


198311 . and it has been 


implemented and tested in some of those popular codes. However, the technique has been 
somew: 

(e. g. 


lat demoted on claims of being notorious 

V slow and less e: 

icient than other methods 

Prsa & Zwitter 

2QQ5; 

Kallrath & Milone 

2QQe; 

Wichmann 

1999 

)• 


Faster methods based on derivatives (LM, DC), start from points near the solution but 
this assumes a previous knowledge of the solution, which is only tuned by the algorithm. 


In this paper we introduce a modified version of the SA method called Adaptive 
Simulated Annealing (ASA), which overcomes the speed problem of standard SA methods 
and ensures fast converging to a global minimum solution for stellar binaries and exoplanet 
radial velocity curves. We have developed an ASA minimization code in IDlQ, which we 
make available to the community, and which can be easily implemented as part of any 
custom binary or exoplanet analysis software. We also provide a detailed guide on how 


^IDL is a high-level commercial programming language and environment by Exelis Visual 
Information Solutions, http: //www. exelisvis. com/ProductsServices/IDL. aspx 
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to implement this ASA method for htting radial velocity curves to a set of radial velocity 
measurements. 

We describe the ASA method in section [2l In section |3] we describe the function that 
the algorithm aims at minimising, and the algorithm itself is described in section H] and [5l 
We present some examples of the performance of the algorithm in section El and discuss the 
results and pros and cons of the algorithm in sections [7] and [HI The code with all its docu¬ 
mentation is publicly available at http://www.cefca.es/people/~riglesias/index.html, 


2. (Adaptive) Simulated Annealing 


SA is a function-evaluation-based technique frequently used in electrical engineering, 
image and signal processing, and other helds to address multi-variable minimization 
problems in which the space dimensions and the complexity or non-linearity is too high 
for conventional minimization algorithms, for instance, those based in derivatives. SA is a 


generalization of a Monte Carlo method i nitially developed for ex amining th e equat i ons o f 


state and frozen s t ates o f n-b ody systems ([Metropolis et ah 


Kirkpatrick et al. 


(19831) and 


19531) . Later on. [Pincus[ fll970[) : 


Cerny fjl999l ) independently generalized the SA idea to solve 


discrete optimization problems that involved local parameter search procedures. 

The initial idea behind SA consisted on developing an algorithm that would 


calculate the (glo 


([Metropolis et al. 


3al) in ternal energy state of a crystalline structure in metallurgic studies 


1953[). The goal of the original SA algorithms was to minimize an 


objective function which resembles the internal energy of materials. The algorithm took an 
independent parameter called temperature and reduced it following a certain annealing law. 
For each temperature, a number Ng^n of possible test states were generated from a previous 
state, and the objective function was computed for each of those test states; a test state 
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was accepted or rejected based on the acceptance rnles, which depend on the temperatnre 
parameter. If a test state was accepted, it became the cnrrent state and, if it was the one 
with the minimnm absolute energy, it was saved. The temperature was then reduced to 
begin the process from the last accepted state. At each temperature, the loop in which the 
states are generated and tested for acceptance is commonly called the Metropolis loop. 

The key feature of the SA algorithm is the way in which the moves to a new test 
state are accepted. All the test states with a energy (objective function) lower than the 
current energy are accepted, but those test states with a higher energy can be accepted or 
rejected based on a probability computed from the energy difference between the current 
and the test state. This is how SA algorithms avoid stopping at local minima, since a state 
with a higher energy configuration, an uphill movement, can be accepted given a certain 
probability. 


Applied to a general minimization problem, the SA method consists of three functions: 
1) a probability density function for a N-parameter space, where N are the parameters to 
minimize. This function generates new test states, 2) a probability of acceptance function, 
which determines whether a given step solution is accepted or discarded, and 3) the 
annealing schedule, which defines how each parameter changes with each iteration step. 


Several yariations of SA are de s cribed in the literature (see e. g 


Corana et al. 

1987; 

Dreo et al. 

2006: 


Otten fc van Ginneken 


1989 


Ingbei 


Salamon et a l 


996 L 


200211 . 


but here we focus on the original SA and the so called Fast Simulated Annealing (FA, 


Szu and Hartley 


19871) techniques. 


Annealing (BA, 


MetroDolis et al. 

1953 

Kirknatrick et al. 

1983) 


mnerally described as Boltzmann 


19831) . Based on the BA notation. 


the probability of acceptance function can be expressed as: 


P = exp{—AE/T), 


( 1 ) 
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where AE is the energy difference between states and T is the schedule of the annealing, 
which is denoted by T by analogy with the temperatures of the different energy states in 
thermodynamics: 


T{k) 


To 

\n{k )’ 


( 2 ) 


In this equation, Tq is the initial temperature, chosen for the algorithm to explore the 
full range of the parameters to be searched, and k = 0,1, ■■■,N with N large enough to 
achieve a convergence. The BA algorithm picks the range of the parameters to be searched 
and starts iterations following this temperature law in search for the ’global minimum 
energy state’ of the system. The generating function for the new states is a gaussian 
function centered in the current state. 


A variation of the SA algorithm is called Fast Simulated Annealing (FA, 


Szu and Hartley 


19871) in which the gaussian generating distribution is replaced by a Cauchy distribution. 


This lead to a better access to distant states due to the long tail of the Cauchy distribution, 
thus improving the exploration of the parameter space in the search for the global minimum. 
To ensure the convergence properties of the new algorithm a change in the schedule of the 
annealing temperatura is needed, following a faster function T(k) = To/k, from which the 
FA algorithm takes the name. Our code implements a generating distribution similar to 
FA, as described in section 12.1.11 


2.1. Adaptive Simulated Annealing 


The rvfit code that we present in this paper to perform fast fitting of rad ial velo city 


curves is b a .sed on he Adaptive Sirn ulated Annealing (ASA) algorithm. ASA (jlngber 


1989I. 


1993 


19961: 


Chen fc Luk 


19991) . was created with the objective of speeding up the 

















convergence of standard SA methods. 


ASA has been already applied to the com pntation o 


with radial velocities and visnal measurements (jPourbaix 


orbits for binary systems 


19981) . But the cases where 


the only observations avaliable are the radial velocities are many more, e.g. exoplanets 
and single/double line eclipsing binaries. In the majority of these cases there are not 
relative positional measurements. Also in these cases, even if light curves exist, sometimes 
is necessary to £t the radial velocities alone to feed initial sets of parameters to more 
elaborated modeling codes. The advantage of this approach is to reach a full solution for 
the physical parameters of the binary system faster. The rvfit code that we present in 
this paper to provide fast fitting of radial velocity curves is based on the ASA algorithm. 


The basic structure of the ASA algorithm is the same of the classical SA. There 
are, nevertheless, some key differences: new distributions for the acceptance and state 
functions and a new annealing schedule; the use of independent temperature scales for each 
fitted parameter and for the acceptance function; and the use of a re-annealing at specific 
intervals. We explain each one of those steps in the following subsections. 


2.1.1. The state function 

Defining U and L as the vectors with the upper and lower bounds of the parameter 
space, in which each point is represented by a vector x, for each parameter Xi the function 
used to generate new test points is: 


xT" = X, + q,{Ui - Li), 


(3) 


where qi is the following distribution function defined in the interval [-1,1] and centered 


around zero 
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qi = sgn{ Vi - -] Ti^genih) 



- 1 


(4) 


In this function, Vi is a random number in the [0,1] interval and Ti^gen{ki) is the 
generating temperature for the parameter Xi, which depends on ki called annealing times, 
an independent index for each parameter. The function sgn{) determines the sign of the 
expresion inside the parentheses. 


We show the behaviour of this distribution function in Figured] for three different 
generating temperatures. The distribution of points given by this function concentrates 
around the central value when Tj ^en is reduced. But even in the case of very low 
temperatures, some points are generated far away from the central point allowing for a scan 
of the space parameters and moving away from the central value whether or not a better 
solution is found. This avoids local minima. 


2.1.2. The annealing schedules 

For this generating function, the asociated generating temperature follows an 
exponential function: 


Ti,gen{ki) =Ti^gen{^)e , (5) 

where Tjgen(O) is the initial value, ki is a natural number which depends on each 
parameter to £t, n is the number of parameters to fit, and c is a constant which depends 
on each problem and need to be adjusted. A similar equation is needed for the acceptance 
temperature: 
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Fig. 1.— Behaviour of the distribution function qi for three different generating temperatures 
Tgen'- dashed for Tg^n = 1000, dotted for Tg^n = 0.01 and solid for Tg^n = 10“^°. The effect 
of reducing the temperature is to concentrate the generated points around the central value 
gi = 0, but still allowing the generation of some points far away. The parameter c was hxed 
to 20 as in our code. 
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Ta = T„(0)e-'="“ 


1/n 


( 6 ) 


Based on the statistical properties of the algorithm f ingbe: 


19891) . a global minimum 


can be reached with these annealing laws while maintaining a fast convergence rate. The 
ASA algorithm becomes adaptive by the dehnition of an independent temperature value for 
each parameter. 


2.1.3. Re-Annealing 

The last difference between the SA and ASA algorithms is re-annealing. The idea 
behind re-annealing is that the rate of change of the annealing schedule can be changed 
independently for each parameter throughout the convergence proccess. In its way towards 
the global minimum, the algorithm travels the parameter space through points with 
very different local topology, i.e. while some parameters vary rapidly in some regions of 
the parameter space, others may vary very slowly in those regions. By using the same 
generating temperatures (Tj ^en) for all the parameters, it is accepted that the cost function 
behaves isotropically, i.e., the topology of the parameters space is nearly the same in all 
points and in all directions. This situation results in a waste of computational effort. 

Therefore, to optimize computational time, ASA decreases T, along the directions 
in which the sensitivity of the cost function is greater, to allow the algorithm perform small 
steps, and incresases R^gen along the directions with a small sensitivity to allow for large 
jumps. ASA adapts its performance by re-scaling the generating temperatures (Ti^gen) and 
the acceptance temperature {To) every Naccept acceptances. 

To adapt the algorithm’s performance, first the sensitivity of each parameter is updated 
based on the local topology of the parameter space. This is done by numerically computing 
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the derivatives: 


dE (Xfoest) 


E{X}yesl “|- ( 5 ^) 

dxi 




where 5i is a small step size in the parameter Xi (see 
generating temperatures are then computed as: 


Chen fc Luk 


19991). 


(7) 


The new 


Ti^geni^^i) \ new 


-T- (k- 


( 8 ) 


Likewise, the new acceptance temperature is reset to Taika) = E{xbest) and Ta(0) is 
reset to the value of the last accepted cost function value. With the new Tj^gen values the 
corresponding annealing times ki and ka are then re-computed as: 


ki 


-log 

c 


- 1,gen 


( 0 ) 




n 


(9) 


k 


a 



r.(o) 


n 


and the algorithm resumes with the new computed values. 


( 10 ) 


3. The rvfit objective function: Radial Velocity Keplerian Orbits 

The model that the rvf it code tries to £t is the radial velocity keplerian orbits equation 
for the case of a double-lined binary, or a single-lined binary where the second object can 
be either another star or a planet. Assuming gaussian uncertainties in the measurements, 
the function that our implementation of the ASA algorithm aims at minimizing is In 
the most general case of a double-lined binary, the function is given by 
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Ni 


x" = Yl 


2=1 


'^calcA 


(Ti 


Nj 

i=i 


'^calcj '^j 


CTo 


( 11 ) 


where Ni and Nj are the number of radial velocity measurements of the primary and 
secondary components, and n* and Vj, and at and aj are the measured velocities of each 
components and their associated uncertainties. In addition, Vcaic,i and Vcaicj in this equation 
corresponds to the expected radial velocities of the primary and secondary components in 
our physical model calculated using the keplerian orbit equation for each component, i.e. 


Vcaic,i = j + Ki[cos {9 + uj) + ecosa;], ( 12 ) 

and 

Vcaicj = 7 + -ft' 2 [cos {9 + u') + ccosw'], (13) 

where u' = u + tt, since the argument of the periastron of the secondary component 
differs by vr from the argument of the periastron of the primary, 7 is the center of mass 
velocity of the system (systematic velocity) measured from the Sun, 9 is the true anomaly, 
e is the eccentricity of the orbit, and Ki and K 2 are the radial velocity amplitudes for the 
two components. 

In equations [12] and [131 th® argument of the periastron of the system, u, is a constant, 
and the true anomaly, 9, is the parameter that varies with time. 9{t) is given by the 
equation: 


9{t) =2 arctan 


1 + e 
1 — e 


tan 



(14) 


where E is the eccentric anomaly, which is computed from the mean anomaly M by 


resolving the Kepler equation 
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E — esiYiE = M (15) 

The mean anomaly M is obtained directly from the orbital period P and the time of 
periastron passage Tp of the system: 


M = -p (( - T,) (16) 

Thus, for a double-lined binary the set of parameters to minimize for a given dataset 
is [P,Tp,e,u,'j, Ki, K 2 ]. For a single-lined binary, where we can only measure the radial 
velocities of one of the components, the set of parameters to minimize is [P,Tp,e,u,'y, Ki]. 


3.1. Resolving the Kepler equation 


The core step in the computation of the objective function is the iterative calculation 
of equation [15] to obtain the value of the eccentric anomaly E. This is a major step in the 
evaluation of the objective function since the Kepler equation must be solved numerically 
and, if not properly done, the evaluation of the objective function will be slow thus 
degrading the performance of the SA algorithm. Since the SA algorithm explores the full 
range of eccentricity and mean anomaly values for a given binary, the calculation of E must 
be robust and must converge quickly for all values of e and M. 


Several routi nes has 
e.g. the review bv iMeeus 


jeen developed in the past to do this in an efficient fashion (see 


199811 ■ Well-known clas s ical me 


1998 

), and the binary search method ( 

SinnottI 

1985) 


hods to solve the Kepler e quation 


19981) . the Kepler method f Meeus 


In the first versions of the rvfit code we implemented a classical Newton algorithm. 
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but soon realized that this method can be improved and accelerated. We therefore ran 
extensive tests to choose a fast and reliable method to solve the Kepler equation for any 
parameter combination. We hnally chose the following method, which we call Newton 
two-step method: 

• We compute an initial value of E using a third order Maclaurin series expansion of E 
and Kepler’s equation. This ensures that the initial value of E is close to its real value 
and therefore provides faster convergence. The implemented equation is the following: 

2 3 

E = M -\- esinM + — sin2M + — (3sin3M — sinM) (17) 

2 8 

• In a second step, a Newton classical iterative method was used to refine the computed 
initial value to the desired accuracy. The E value previously computed is used to 
feed this iterative method. At each iteration a new value of the eccentric anomaly is 
computed using the equation: 


r-, r-, M + e sin Efc — Efc 

^k+1 — -Wk -\ - - - p; - 

1 — e cos Ek 

This process is repeated until \Ek+i — Ek\ < e. In this case the accuracy was set to a 
somewhat high value of e = 10“^*^ radians to avoid systematic errors. 

The number of iterations versus the eccentricity of our Newton two-step routine and 
other classical routines is shown in Figure |2] for an M=0.3 radians value of the mean 
anomaly. This value was deliberately selected to show some strange effects seen in the 
number of iterations at high eccentricity values for some classical routines. The discrete 
jumps in the plots are due to the fact that iteration numbers are integers, so jumps appears 
for small values in the logaritmic vertical axis. The comparison among our two-step method 
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and the other methods are displayed in the left panel of the Fignre[2]for the fnll range of 
eccentricities. 


Since we selected 32 steps in the binary search method, the plot of the iteration nnmber 
is a constant. In the Kepler method we also implemented the two-step approach mentioned, 
computing the initial value of the eccentric anomaly with the eq. [T71 The Kepler method 
needs from 3 to 30 iterations depending on the eccentricity value but the behaviour shows 
a regular trend with a maximun of iterations near e=0.7. The Newton classical method is 
very fast, converging to the solution in less than 10 iterations for eccentricities below e=0.95 
but shows a strong increment near the limit of e < 1. To investigate this behaviour we run 
aditional tests for the eccentricity in the range [0.95,0.999] for each method. The results are 
displayed in Figure |2l A strong random variation is observed in the number of iterations 
needed to achieve the desired accuracy, with peaks above 200 iterations for some values 
of the eccentricity and less than 20 iterations for the neightbour test. No clear pattern is 
observed in this variation, which suggest a weak convergence in this range of eccentricities. 
The results of the tests shown here are for M=0.3, but we hnd a similar behavior for other 
values of M, affecting slightly different intervals of eccentricity near the limit of 1. 


has observed this same effect (see his Figures 4 and 5) and claims that 
the problem is present for values of M near 0 and eccentricity near 1, i.e. the parameters 
that makes the denominator in the equation [TSl almost 0. He quoted some solutions, which 
are reduced to limiting the value of the eccentric anomaly correction in the algorithm. 


MeeusI (19981) 


Our implemented Newton two-step algorithm avoids this problem and improves the 
performance of the Newton method by starting from a solution very near the correct result, 
thus reducing the number of iterations and increasing the stability. 
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4. The implemented ASA algorithm 


In this section we provide a detailed description of our ASA al gorithm implementa tion 


For our implementation we follow the same approach detailed by 


Chen &: LukI ( 19991) . 


except for how we calculate the initial acceptance temperature and other improvements, as 
will be described in more detail below. The code was implemented in IDL 7.0. The ASA 
code is packed in a main program called rvfit but several common routines are separated 
in other file for convenience. 


Based on the equations described in the previous section, and assuming no previous 
knowledge of the orbital parameters of the system, our ASA algorithm minimizes the 
function in eq. [TT] in the multidimensional space given by the following set of parameters: 
the orbital period of the system, P, the time of periastron passage, Tp, the eccentricity of 
the orbit, e, the argument of the periastron, u, the systematic velocity of the system, 7 , 
and the radial velocity amplitudes, Ki and K 2 (only Ki in the case of single-lined binaries 
and exoplanets). 

ASA works by assuming that each measurement is in the form of [ti, RVi, a{RVi)], 
i = 1, ■■jNi, where ti is the time of observation i (in HJD or BJD), RVi is the measured 
radial velocity (m/s or km/s), and a{RVi) is the uncertainty of that radial velocity in the 
same units as RV. If the binary is double-lined then a set of parameters [tj, RVj, a{RVj)], 
j = 1,.., Aj, for the secondary need to be provided. We assume those uncertainties follow a 
gaussian distribution. 


4.1. Initialization 


Our algorithm implements some improvements over the one by 


Chen fc LnkI ( 19991 ) 


One of those improvements is the initialization of the acceptance temperature {Tao). 








we use the 


Chen fc Lukl (1199911 set T„n fr om the value of the cost function in an aleatory state. Instead, 


Dreo et ah 


(120061) description with the following steps: 1) we compute a set of 


aleatory perturbations is computed, i.e. 100, 2) we compute the average of the cost function 
variations due to these perturbations < >=< {xl+i ~ xi) 3) we stablish an initial 


acceptance probabilit y oi v = 0.25, claime c 


parameter spaces (see 


Gelman et ah 


2003 


to be 


Driscoll 


le optimal value for 5-100 dimensional 


2006 


MacKay 


20061), and hnally 4) we 


compute the optimum Tao from p and < > as: 

1 


P = 


1 + ^<Ax^>ITaO 


TaO = 


< Ax^> 
ln(l/P + l) 


(19) 


This change ensures an appropriate initial acceptance temperature Tao, optimizing 
convergence and therefore reducing computational time. 


4.2. Parameter tunning 

An advantage of the ASA algorithm is that some of the control parameters are 
automatically and easily set during the first stages of the execution and are tuned along 
the evolution of the temperatures. However, a few of them, namely Nacep, Ng^n and c, need 
to be set by hand and adapted to each particular problem.. Here we explain how we adapt 
those parameters in rvfit. 

The c constant (introduced in eqs. [5] and [6]) controls the annealing rate, at which the 
generating and acceptance temperatures are reduced. When c is too low the algorithm 
needs more time to achieve a low temperature stage where the neightbourhood of the global 
minimum is explored and refined. So the algorithm slows down. When c is too high the 
temperature drops faster and the algorithm can get stuck near a local minimum depending 
on the cost function and the parameter space topology. 


Chen fc Lukl (119991) claim values of c in the range [1, 10] work well. However, our tests 
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with published observations conclude that c = 20 provides faster convergence and works 
well for the number of parameters we need to £t for. This large value of c makes mandatory 
the use of double precision arithmetic in the code variables to avoid computational overflows 
or underflows. For some cases, it may be convenient to reduce this value. 


We set the other two paramet ers that need t unnin g to W^en = lO'^ and Naccept = 10^, 


an order of magnitude higher than 


Chen &: Luk fjl999l) . These values are justihed by the 


need to have better sampling of the orbital period in cases where the period of the system 
is small compared to the range of periods the code explores. A good example is when short 
period binaries are found in long observing campaigns. Of course, these values can be 
revised for some cases, but we find they work well for a wide range of synthetic and real 
radial velocity datasets (see section ED. 


4.3. Termination 


Another change implemented in our rvfit code is the way in which a termination 
condition is ac hieved. A number of stopping rules are described in the literature (see 


Locatelli 


20021. p. 23), all of them based on the fact that when the algorithm does not evolve 


signihcantly over a number of temperature changes then it must be stopped. When fitting a 
model to a set of observations, some aditional information is provided by the uncertainties, 
cr: assuming that the model is correct and the uncertainties are well computed, then the 
fitted model must pass through most of the points at a distance less than cr. When the 
algorithm is stuck near a local minima we can use this criterion wheter it is stopped. 


We decided to use a mixed criterion based on the information provided by the 
uncertainties and the progress of the algorithm. Our criterion gathers information about 
the convergence state as the algorithm evolves, and with that information and the actual 
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value of the makes a decission about stopping. To do this, we make a list with the last 
best values of (in our case iVg = 5), and append the best fo the list after each 
temperature step. When the differences among all the values of iu the list are less than 

£ (with e: = 10“^), the annealing loop is stopped if the value of the last x^ is: 


x" < E 


k=l 



( 20 ) 


where a = (X]r=i is the mean of the uncertainties ai in the data (primary and 
secondary), with n = Ni + Nj. This upper limit is set assuming that in a good fit some of 
the measurements are separarated less than a from the model. 


If this last condition is not met, we reset the acceptance temperature = T^o, and the 
program resumes the annealing loop. 


To prevent the loop from running indefinitely, we impose a limit, Werm, the number of 
reannealings that the annealing loop can do without an improvement in the x^ value. This 
condition would only happen if the algorithm got stuck in weird values, when the initial 
values provided makes the algorithm unstable, or if the htted data doesn’t correspond to 
the model. We impose a high value of Werm = 20 to allow the code to achieve a stable state 
before termination. 


After terminating the annealing loop, the code computes the uncertainties in the 
parameters and uses the parameter values to derive the physical quantities of the system. 

The behaviour of the algorithm is shown in Figures [3] and HI Those hgures show 
the acceptance temperature and the cost function values for all the accepted parameter 


conhgur ations for a fi 


LV Her flTorres et ah 


ting 


run to the data of the eccentric double-line eclipsing binary 


20091) . In this fit all the parameters were leave free, including the 


orbital period. The figures illustrate the evolution of the acceptance temperature with the 
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annealing phases of decreasing T^, the snbseqnent reannealing jumps, and the associated 
groups of new acceptances with higher cost function values when a reannealing occurs. 
Figure O shows the resulting radial velocity curve £t for this binary. 

Since ASA is a heuristic method with aleatory generation of test parameter values, 
these plots will be somewhat different each time the code is executed. However, they will 
have the same overall behaviour of Tq and due to the annealing-reannealing cycle. 


5. Computation of uncertainties 

Once the algorithm has converged to a global minimum solution, rvfit computes the 
uncertainties of each parameter in one of the two ways described below. 


5.1. Fisher matrix 


The Fisher matrix is a popular procedure to obtain quick uncertainties in the model 
parameters. It is based on the central-limit theorem which states that ’a well-behaved 
likelihood distribution is asymptotically Gaussian near its maximum’. As a consecuence 
this procedure can only describe elliptical uncertainty contours. Their implementation is 
straightforward and numerical derivatives can be used if the function doesn’t have analytical 
form. 


Following 


Andrad (1201011 . rvfit computes the Fischer matrix of all parameters at the 


global minimum position as: 


IF] 


c 52 

92 

dx‘‘ 

dxdy 


92 

dydx 

9 j /2 




X 
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The computation of the partial derivatives is done following 


Coel ( 2 OO 9 I). Once the [F] 


matrix is computed, the covariance matrix is obtained by inversion: 


^xy 1 

( 22 ) 

'^yx ) 

The main drawback of computing the uncertainties with the Fisher matrix is that 
it cannot represent non-linear correlations between parameters, which can be seen in a 
number of models. 



5.2. Markov Chain Monte Carlo uncertainties 


rvf it also includes a Markov Chain monte Carlo (MCMC ) method for the cases 


where more detailed uncertainty distributions P{x) are desired. 


Ford! (12005 


2006 ) applied 


a bayesian analysis to the computation of uncertainties of exoplanet radial velocity curves. 
Much of his work can be directly applied to this problem since the Metropolis-Hastings 
(MH) algorithm he used to sample the P{x) distribution uncertainties is the same. 

In our case, a separate MCMC code to compute th e uncertainties was implemented. 


using the simple MH algorithm depicted in 


Fordl (120051) . This code works in two basic 


steps. First, a new sample x' is generated from a candidate distribution q{x). In our case, 
this distribution q{x) is a multidimensional gaussian function of width a, which is a vector 
composed of the widths for the distributions of each of the htted parameters. Second, 
a decision is made to accept or to reject the new state, using the Metropolis-Hastings 
acceptance probability A, given by 


/ P{x') q{x,x')\ 

V ’ P{x) q{x', x) J 


A(x, x') = min 


(23) 
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In the cases where the distribution q{x) is symmetric the values q{x,x') and q{x',x) 
turn out to be equal and the MH algorithm is the called Metropolis algorithm, the same 
one used for the SA algorithm. This is the situation in the typical case where a gaussian 
distribution of width a centered in the actual state x is chosen, as in our MCMC code. 
Based in the ratio of likelihoods of the two states, the new acceptance probability is then 
given by: 


A{x,x') = min ( l,exp 


;(,P(I)-XV)) 


(24) 


The Markov Chain is built with the accepted values x' , if these were accepted, or with 
the old X, if there is not an acceptance. The old x must be appended to the chain whenever 
there is not an acceptance. If this is not done, some bia s will be introduced in parameters 


with a hard bound like e = 0. See flEastman et ah 


20131. sec 3.5.3) for a discussion. 


Our MCMC does not make use of a burnin phase prior to the computation of the 
Markov Chain because the chain is started in the parameter values provided by the ASA 
algorithm which are very near the maximum of the unknown distribution. But the width a 
of the proposal distribution is a key parameter which needs to be tuned up to achieve a good 
mixing of the Markov Chain and to obtain a fast sampling of the objective distribution. 
This a is usually selected depending of the desired acceptance probability to achieve a 
good mixing and we selected again an acceptance probability of p ~ 0.25, as in the TaO 
initialization (see section ITTl) . Given that each htted parameter has its own range and 
marginalized distribution the a values must be unique for that parameter. 


Our procedure to tune up this vector is as follows. To tune up this vector, we hrst set 
the initial a to 1/3 of each parameter range, and we measure the acceptance rate for 10^ 
generated states by computing the individual acceptances using eq. [2ll If the acceptance 
rate is less than p ~ 0.25, a is divided by 1.2 and the process begins again in the previous 
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point. This process is repeated until the acceptance rate exceeds the pro posed 


value. We are aware that other more elaborated algorithms exist (see e.g. 


Ford 


0 ^ 0.25 


2006 


sec 


3.2), however this simple algorithm is fast and does a good job providing a good mixing of 
the Markov chain. 


Once the a for the proposed distribution is computed, the MCMC algorithm computes 
the uncertainties for all the parameters. Usually, a MCMC chain with 10^ points is enough 
to provide a reliable measurement of the uncertainties, but longer chains might be needed 
in some cases. 


6. rvfit perfomance tests 


To test the performance of rvfit we compared its results to published solutions for a 
set of systems. The systems were selected to cover a range of physical and observational 
conhgurations, i.e. double-lined and single-lined binaries, with both circular and eccentric 
orbits and with long and short periods, and exoplanet systems. We also selected datasets 
with a variety of observing conditions: with few and many observations, with high and low 
signal to noise measurements, and densely and loosely observed. 


We do not include datasets aimed at observing spec ific effects, such as the Rossiter- 


McLaughlin effect 


Svbilski et ah 


201811 . 


see 


Ohta et ah 


2005 


Gimenez 


20061 ). or relativistic or tidal effects (see 


We selected two double-line spectroscopic and eclipsing binaries, two single-line 
spectroscopic binary stars, and two exoplanet systems. The published solutions for those 
systems obtained from the literature and the results of our rvfit code are shown in Table 
[U All uncertainties were computed using a MCMC chain with n = 10® samples, except for 
GU Boo with all parameters free (n = 10® samples) and HD 37605 (n = 5 x 10® samples). 
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6.1. LV Her 


LV Her is a double-line eccentric eclipsing binary star flTorres et al.l 120091) . We selected 


this system because its high eccentricity, the high quality of its radial velocity curve, the 
long orbital period, and the long observing time interval relative to the period. For this 
system two fits were done: the first one leaving all the para meters free and the s econd one 


fixing P to the value from the light curve to mimic the fit by 


Torres et ah 


(120091). 


Our first fit arrives to a solution remarkably similar to the published one (see Figure 
|5]). Our fitted period has a greater uncertainty since the one published was derived from the 
photometric eclipses of the system, which are high quality time marks. Even so, our fitted 
parameter set, based only in the radial velocity data, recovers perfectly the overall solution. 


Fixing the orbital period to the value from the light curves ((Torres et ah 
again, similar parameters to the ones published. 


200911 yields. 


6.2. GU Boo 


GU Boo is a double-line and eclipsing binary star (jLopez-Morales fc Ribas 


20051), 


composed of two nearly equal low-mass main sequence stars (Mi _2 — O. 6 M 0 ) in an circular 
orbit with a short orbital period of 11.7 hours. 


Leaving all the parameters free the eccentricity obtained with rvfit is compatible 
with zero and the uncertainty in u is fairly large. This suggest a circular orbit. The MCMC 
run shows fairly large uncertainties in those parameters related with the periastron: the 
marginalized histograms for Tp, e, and u are severely spread over a wide range and those 
of Tp and u are multimodal and runing a longer MCMC chain doesn’t fix the situation. 
This behaviour is a consequence of the degeneracy in Tp and u due to the null eccentricity 
and is a excellent indicator of this situation. In Table [1] we quote for P the uncertainty 
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derived from a gausian fit to the histogram and for the other parameters the uncertainties 
and values derived from the maximum of the histogram and the 68.3% shortest conhdence 
interval. We show the htted RV curve and the residuals in Figure E] 


In the second fit we fixed P, Tp, e and oj to the values derived by 


Lopez-Morales fc Ribas 


(120051) from the system’s light curves, u was fixed to 90 degrees to match the time of 


periastron with the instant of the primary eclipse. This leaves only three parameters 
free which were easily recovered in our £t. The uncertainties are gaussian-like and larger 
than the published values since in the original paper they were computed by simple error 
progagation. 


6.3. GJ 1046 


GJ 104.6 is a sy stem containing a brown dwarf in an eccentric orbit around an M dwarf 


flKlirster et ah 


20081) . The system is single-lined because of the low luminosity of the brown 


dwarf secondary. We fitted all the parameters simultaneously and obtained remarkable 
agreement with the published values (see Figure [7]), obtaining a smaller value for the 
than the published value. The MCMC uncertainties are very well modeled by gaussians 
and in Table [U we quoted as uncertainties the a of those gaussians. 


6.4. HR 3725 


HR 3725 is a sing le-lined spectroscopic binary containing a G2 III primary 


flBeavers &: Salzer 


19851) . Presumably, the secondary is a dwarf star with a late spectral 


type in a circular orbit. Using the published dataset we htted two models: the hrst one 
with all the parameters free, which confirms the hipothesis of a circular orbit. In the second 
model we hxed e = 0 and oj = 90, as the hrst model suggest. In both cases, our period is 
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smaller t 
stated in 


ran the publishec 


Beavers fc Salzer 


value , which was obtained by an independent algorithm not 


fjl985[) . The Ki value agrees within uncertainties. Our fitted RV 


curve is shown in Figure [Hi 


6.5. HD 37605 


HD 37605 was the first exoplanet discovered by the Hobby-Eberly Telescope 


flCochran et al.ll2004i) . The planet has a very eccentric orbit (e = 0.737) with a period 


of 54 days and a mass of msini = 2.84 Mj. Again, we did two fits, one leaving all the 
parameters free, and another fixing 7 = 0 to mimic the published fit. For both models, the 
resulting parameters are in good agreement with the published orbital solution. The htted 
RV solution is shown in Figure [9l 

For this system, the MCMC was run over 5 x 10^ times, due to the asymmetry in the 
marginalized histograms for the model with all the parameters free. 


6.6. Kepler 78b 


Kepler 78b is the first Earth-like transiting exoplanet discoyered in the Kepler mission 


data flPepe et al. 


2013 


Sanchis-Oieda et al. 


2013 


Howard et al. 


20131) . It orbits its parent 


solar-like star in a nearly circular orbit and, due to their small mass-ratio, the HARPS-N 


dataset (jPepe et al. 


20131 ) has a very low signal to noise ratio (see Figure fTOj) . So, this 


system poses a very stringent test to our code since it must recover the very low-amplitude 
RV signal caused by this exoplanet. 


In a first £t, we kept Tp, e, u fixed following a similar approach to 


Pepe et ahl (12013! 


Tp was set to the reference epoch in their Table 1, the eccentricity was set to 0 and oj was 
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set to 90 degrees to match the instant of the transit with Tp. Note that they htted the mean 
longitnde Aq at epoch Tq. Onr £t results in slightly lower values of P and Ki, near the limit 
of uncertainties. In a hrst test we set a period range between 0.1 and 1.0 days, but the 
algorithm found solutions with P = 0.265024 (x^=327.73) and P = 0.748368 (x^=330.62), 
both with lower than the Kepler period (y^=339.98). The presence of those solutions 
led us to constrain the period search to the interval [0.3, 0.5] days to match the period of 
the planet from Kepler data, since this is a local minimum in a wider period domain. 


Since there is a precise period computed by 


Sanchis-Oieda et af 


( 20131) from the Kepler 


mission photometry we ran another £t but no w fixing P, e, anc 


better agreement with the value published in 


Pepe et ah 


u. The value for Ki has a 


( 2013 ). This is due to the small 


change in the period which varies the overall solution, as shown in the figure. 


7. Discussion 

7.1. Treatment of circular orbits 

As shown in section [6] when the eccentricity is null the periastron remains undehned, 
so the parameters related can be set to a convenient value. In eclipsing binary stars or 
transiting exoplanets with circular orbits it is common to dehne the periastron as the point 
of the inferior conjunction since it is well constrained by the center of the main eclipse (Tq) 
and can be used as the phase reference. So it is common to £x Tp = Tq and u = ti/2 to 
match the dehned line of nodes with the vision line. 

This situation is easy to identify when the code is forced to £t the eccentricity in a 
circular orbit, since the uncertainties in the fitted values of Tp and u: are quite large, and 
the htted e value is near 0 with a large uncertainty. In this case, the best procedure is to 
re-£t hxing e = 0. 













Object Name 

P 

Tp 

e 

UJ 

7 

Ki 

K2 

X3 

Remarks 

Reference 


(d) 

(HJD or BJD) 

- 

(deg) 

(km/s) 

(km/s) 

(km/s) 




Double-line binary systems 

LV Her 

18.4359535'f 

2453652.19147'^ 

0.61273(73) 

352.20(.24) 

-10.278(94) 

67.24(.19) 

68.59(.27) 

- 

- 

[1] 


18.43600(14) 

2448414.904(29) 

0.6137(19) 

352.22(.27) 

-10.291(88) 

67.31(.24) 

68.67(.33) 

81.416230 

all free 

[this work] 


18.4359535'f 

2448414.9117(61) 

0.6138(20) 

352.23(.28) 

-10.288(94) 

67.32(.26) 

68.68(.35) 

81.498606 

P fixed 

[this work] 

GU Boo 

0.4887280'f 

2452723.9811^ 

O.O-f 

- 

-24.57(.36) 

142.65(.66) 

145.08(.73) 

- 

- 

[2] 


0.488877(73) 

2452724.033(+(®) 

o.oii(j:°®) 

137(1 

-24.81(1;™) 

142.5(11“) 

144.8(111) 

46.712033 

all free 

[this work] 


0.4887280'f 

2452723.gSll'f 

o.of 

907 

-24.64(.69) 

142.4(1.3) 

144.9(1.5) 

51.254343 

P, Tp, e, uj fixed 

[this work] 





Single- 

•line binary systems 





GJ 1046 

168.848(30) 

2453225.78(.32) 

0.2792(15) 

92.70(.50) 

- 

1.8307(22) 

not aplicable 

12.7 

- 

[3] 


168.845(21) 

2453225.78(.30) 

0.2792(12) 

92.69(.49) 

-0.0672(21) 

1.8307(20) 

not aplicable 

12.466580 

all free 

[this work] 

HR 3725 

66.717(4)-f 

2444667.86(.25) 

0.0^ 

- 

-18.93(.38) 

23.30(.52) 

not aplicable 

- 

- 

[4] 


66.6683(+g^) 

2444665(1^®) 

0 .002(j:J») 

344(llf) 

-18.27(1;37) 

23.14(1;«) 

not aplicable 

45.036888 

all free 

[this work] 


66.6673(47) 

2444617.93(.16) 

0.0^ 

9o7 

-18.28(.26) 

23.08(.47) 

not aplicable 

45.056023 

e, tj fixed 

[this work] 

Exoplanet systems 

HD 37605 

54.23(.23) 

2452994.27(.45) 

0.737(10) 

211.6(1.7) 

- 

0.2629(55) 

not aplicable 

- 

- 

[5] 


54.30(j:;32) 

2452939.82(+;^^) 

0.7351(87) 

211.1(1.4) 

0.0002(16) 

0.2618(1“*) 

not aplicable 

42.753639 

all free 

[this work] 


54.21(.17) 

2452939.96(.51) 

0.7342(78) 

211 .0(1.1) 

0.07 

0.2623(49) 

not aplicable 

42.326212 

7 fixed 

[this work] 

Kepler 78b 

0.3550(4) 

2456465.076392^ 

O.O'^ 

- 

-3.5084(8) 

0.00196(32) 

not aplicable 

- 

- 

[6]. [7). [8] 


0.34951(15) 

2456465.076392^ 

0.0^ 

9o7 

-3.51035(21) 

0.00173(11*) 

not aplicable 

339.98215 

Tp, e, uj fixed 

[this work] 


0.35500744(6)'f 

2456464.1575(98) 

o.of 

9o7 

-3.51044(21) 

0.00186(28) 

not aplicable 

326.59471 

P, e, Lo fixed 

[this work] 


Table 1: Results for the fitting of radial velocities data for the selected systems. All the uncertainties were computed using a MCMC chain 
(see text). For each parameter, the symmetric uncertainties were computed by fitting a gaussian function to the marginalized histogram. 
In this case the quoted uncertainty is the a of the fitted gaussian function. The asymmetric uncertainties were computed using the 68.3% 
shortest confidence interval over the histogram. In all the cases, the uncertainties are expressed as the two last digits of the parameter, 
with the decimal p oint present where 


refere nces are: [ 1 ] = 


(1201311 . [8] = 


Torres et al 


(2QQ5), [3]=Kyraigrji^ 


Howard et al 


(2008D, [4]= 


( 201 ^ 


( 200911 fthere is a typ o in th e pu blished value of 7 [To rres, priv. comm.)), [2] = 


priv. comm.!), [21 = 

Louez-Morales & Ribas 

Sanchis-Oieda et al 

( 

201.1 

), [7]= 

Pene et al. 
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7.2. Speed tests and dependence with the number of observations 


A primordial mot ivation for all fitting code s is the execution speed. Some authors 


claim that SA is slow. 


Milone &: KallrathI fj2008f) publi shed an exten d ed an alysis of classical 


SA and its performance using an algorithm similar to 


Cor ana et ah 


( 19871) . and applied it 


to the computation of eclipsing binary stars light curves. But they only fitted light-curves, 
not RV curves, and used the WD code to compute the cost function. Fitting three or four 
parameters, their measured execution times are as long as 6.2 x s, though in one fit they 


get a more reasonable value of 386 s. 


Prsa &: Zwitterl (120051) claim that SA methods ’are 


notoriously slow’, without going into details. In their concl usions they s tate that ASA, with 


Powell’s direction method, is a ’very promising candidate’ 


Prsa 


merits applied to eclipsing binaries ’are still to be confirmed’. 


(2011) claims t 


rat ASA 


Kallrath fc Milonel ( 2006 ). in 


section 4.3.4.3, states for SA that ’the computational cost to reac 


as temperature annealing has to be sufficiently slow’. 


1 this s tate can be high. 


Wiclunann 


( 19991) includes SA in 


Nightfall, and in their User Manual states that ’Be prepared for a computing time on the 
order of a day or more’. We presume that he implemented a standard SA though no details 
are given. 


To our knowledge, these claims were made in the context of fitting simultaneous 
photometric and RV observations of eclipsing binaries where a complex model need to be 
fitted to account for the variety of phenomena in such systems. This triggers the number of 
parameters to be fitted simultaneously, including tidal distorsions, limb darkening, gravity 
brightening, mutual heating effects, computation of visible photosphere due to eclipses, and 
stellar spots. It is clear that in these models the cost function is harder to compute and the 
space parameter is harder to search. In addition, some of these codes implement classical 
SA, not the newer versions with better convergence properties. There have been, however, 
developments in computational technology in recent years, which allow computations with 
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much higher speeds. 

The computation time of the ASA algorithm depends on a number of factors: 

• the number of parameters to £t, 

• the complexity of the functional dependence of the cost function with the parameters, 

• the size of the dataset, and 

• the particular conhguration of internal variables to run the code, particularly c and 
N 

' gen • 

In addition, the ASA is a heuristic algorithm with a strong random behaviour. This means 
that each time it is run the trajectory in the parameter space will be different, in spite of 
starting the algorithm with the same initial parameters. 

Because of those factors, the precise execution time of the ASA algorithm cannot be 
measured or predicted since each run of the code will lead to different values even with the 
same dataset and starting points. Thus, we decided to estimate the execution time of the 
algorithm taking a statistical approach. For a model radial velocity dataset, we ran the code 
a large number of times to obtain statistics on the times needed to complete the £t. These 
speed tests were done to provide a complete physical solution, not only the solution of the 
ASA algorithm. This includes the code needed to read the data hies, to compute the initial 
parameter values and, after arriving to a solution, to compute the physical parameters of 
the system with their uncertainties, and to display the plots with the fitted radial velocity 
curve, the observed data and the residuals. The tests do not include the computation of 
the uncertainties using the MCMC code, which is in separated routines. The speed tests 
were done in a DELL inspiron laptop with an Intel Core 15 processor and 8 GB of RAM 
memory, running Ubuntu Linux 12.04 LTS. 
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To assess the dependence of the execution time with the number of RV observations we 
simulated diferent datasets with a particular orbital conhguration. We chose an orbit with 
the parameters P = 10 days, Tp = 2450000.0, e = 0.1, ca = 90 degrees, 7 = 0 km/s, Ki = 20 
km/s. The chosen orbit depicts a single line system without loss of generality, since for a 
double-lined system the total observed points are distributed between the two stars with 
the same computational load for computing the cost function as in a single-lined system 
with the same total number of points. 

With this conhguration we simulated four datasets with 15, 50, 100 and 1000 randomly 
spaced data points over a span of 3 periods. Gaussian distributed undertainties with 
O'noise = 2 km/s Were added to each dataset, resulting in radial velocity curves with a 
signal-to-noise ratio (SNR) 10. Each model was htted 1000 times using rvf it. For each run, 
the execution time and the was logged. Figure fTTl shows the execution time histograms 
of the four datasets. To quantify this plot we dehned tes; ^ 95 ; and fgg as the execution times 
encompassing the 68%, 95% and 99% of all runs. The results are summarized in Table |2j 
Typical execution times are below 20 s for common single-lined binary datasets and below 
30 s for double-line binary datasets. Earth-like exoplanets with datasets of hundreds of 
points could be analized in less than two or three minutes. 


7.3. Robustness tests 

To make the code as robust as possible, we performed thousands of htting tests by 
generating synthetic datasets with a wide range of orbital parameters varying P, Tp, 7,^1, 
and K 2 in a aleatory fashion and covering the full range of values for e and u. A number 
of flaws were detected, mainly related with overflows and zero denominators, which were 
corrected checking the conditions of the computations in volved. 
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Number of points 

Execution times (s) 


^68 

^95 

^99 

15 

4.1 

6.8 

8.9 

50 

7.6 

12.7 

16.1 

100 

20.2 

30.4 

38.5 

1000 

81.1 

132.7 

175.1 


Table 2: Execution times measured over 1000 consecutive runs for the same synthetic datasets in Figure 
m Each time was measured from the histogram computing the cumulative distribution function (CDF) and 
interpolating the desired percentage. Note that doubling the tes execution times is a good guess for tgg. 


7.4. Refinement of the result 


In th e ir wor k about the ASA algorithm applied to combined visual and RV observations, 


Pourbaixl fjl998h claims that a local search has to be used to tune the minimum found. 
This could be due to the fact that the termination condition for his ASA algorithm stops 
the computation after a fixed number of temperature reductions with the goal in mind 
of obtain a automatic value for the annealing-rate parameter c. Our criterion to stop the 
annealing loop is somewhat different since we keep a track of the variations in the cost 
function to have control over the convergence. So the refinement is an unnecesary step 
for our algorithm. Whether the result is to be refined, the output parameters from the ASA 
algorithm have to be th e input parameters for a code such as Levemberg-Marquadt (see 


Wright fc Howard 


20091). 
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8. Conclusions 

The Adaptive Simulated Annealing Algorithm (ASA) provide a valuable tool to £t 
functions in highly- dimensional parameter spaces. The first versions of the Simnlated 
Annealing algorithm, the so called Boltzmann Annealing, was compntationally slow. The 
new developments in this algorithm, namely Adaptive Simulated Annealing (ASA), makes 
it an option to take into acconnt. 

With present domestic compntational technologies, a complete solntion for radial- 
velocity curves with preliminary uncertainties can be obtained in times of the order of tens 
of seconds or less. In eclipsing binary systems and transiting exoplanets, where the period 
can be determined from their light curves, the fitting is even faster obtaining a complete 
solntion in seconds. 

ASA allows to fit radial velocity cnrves leaving free all the parameters, inclnding 
the period. This is an advantage over others algorithms where the period must be fixed 
beforehand using other techniqnes snch as periodograms. Fast Fourier Transforms or Phase 
Dispersion Minimization techniqnes. 

One of the advantages of the ASA approach is that no derivative is needed to compnte 
the £t since it is based only on fnnction evalnations. This efficiently avoids local minima. 
Also, this approach allows to concentrate all the physics in one fnnction, the objective 
fnnction, in our case the where the physical model comes into the Vcaic- More elaborated 
models can easily be intodnced changing this fnnction. 

Due to these refinements, the advantages of ASA over the classical SA algorithm are 
clearly stated: 1) Better convergence properties in high dimensionality spaces with great 
topological complexity, 2) Individnal adaptation of each parameter to the local topology, 
and 3) Faster schednle for the generating and acceptance temperatnres. 
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But this advantages have some cost. The implementation of the ASA algorithm 
is more complex than the classical SA algorithm, there is a larger number of internal 
parameters, namely, the temperatures and the sensibilities and, 3) the exponential form of 
the temperature schedules force to check the exponents to avoid overflows and force the use 
of double precission aritmethic. 

We have developed a fitting code called rvf it which makes use of ASA, to fit keplerian 
radial velocity curves. Our code, implemented in IDL, computes initial uncertainties using 
a Fisher matrix but we have also tested a routine to obtain uncertainties using a Markov 
Chain Monte-Carlo (MCMC) technique with the Metropolis-Bastings sampler. 

The rvf it code, which is publically avaliable at http: //www. cef ca. es/people/~riglesias/index .h 
shows their full capabilities when a search in the full parameter space is needed. Thus, this 
code may be helpful in situations where the orbital period is unknown, as may be the case 
of RV surveys, since it avoids the previous computation of a periodogram or the use of 
other routine to find the period. Also, it provides the full set of parameters in a consistent 
way taking into account the observational uncertainties of the data. 

Other, simpler techniques, can be applied in the cases where only a small subset of 
parameters need to be fitted. For instance, for a circular orbit of known period, a sine 
curve can be fitted directly to the radial velocities. This is because there is a limit to the 
minimum number of computations that must be done at each temperature to reach the 
termination condition. This mechanism is of low efficiency for a small parameters set, where 
other algorithms can perform better. 

Our tests with real and synthetic radial velocity curves are very promising, with our 
computed values and uncertainties in full agreement with the published values, thus proving 
the power of this technique. 
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We note that the current version of rvfit is only suited for htting single exoplanet 
systems. Simultaneous multi-planet hts will be implemented in future versions of the code 
Future rehnements of rvfit will also introduce Keplerian perturbations in the orbit due 
to the presence of other bodies in the system, such as interacting planets or a thirth star 
in an external orbit, calculations of jitter contributions for exoplanet RV datasets and the 
capability to merge and simultaneously fit RV datasets from different sources. 
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Fig. 2.— Number of iterations versus eccentricity for four different methods to solve Kepler’s 
equation. The left-side panel shows 100 evenly distributed eccentricity values between e = 
0 — 0.9999, for M=0.3 radians. This plot reveals how the classical Newton method has 
convergence problems when e ~ 1. The right-side panel shows the eccentricity interval 
e = 0.95 — 0.9999, where we compute 100 values for M = 0.3 radians. Both plots show how 
our implemented Newton two-step method is more efficient than all the others. 
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Fig. 3.— The typical behaviour of the acceptance temperature Ta (upper panel) consists of a 
sucession of annealings and reannealings. Each upward jump in the temperature corresponds 
to a reannealing and it is associated to the exploration of an extended region in the parameter 
space, which turns in new acceptances. Those acceptances can be seeing in the plot of the 
cost function (lower panel) as upward spikes, since the new acceptances typically have a 
higher cost function. This behaviour allows the ASA algorithm to avoid local minima. 
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Fig. 4.— Behaviour of the parameter values during a run for the LV Her dataset. In the 
amplitudes plot, the black line represents the acceptances for Ki and the red one for K 2 . Note 
the high density sampling of the parameter domain in the hrst stages associated to a high 
temperatures phase (both and Tj ^en) in the algorithm. Once the global minimum region 
is found, near the acceptance 3400, all the acceptances are produced in that neightbourhood 
and the parameters evolution is estabilized. Also, note that the hrst parameter to become 
stable is the period, which show huctuations smaller than the thickness of the line in the 
plot. 
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Fig. 5.— Radial velocity curve for LV Her for the fit with all the parameters free, including 
the period. The error bars are smaller than the point size in the RV plot. This solution is 
visually indistinguishable from the one with P hxed. 
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Fig. 6.— Radial velocity curve for the low-mass eclipsing binary GU Boo. This plot corre¬ 
spond to the second £t, with a circular orbit and in which we fixed P and Tp to the values 
obtained from the photometry. The fit with all the parameters free is indistinguishable from 
this plot. 
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Fig. 7.— Fitted radial velocity curve for the single-lined binary GJ 1046. For this plot all 


the six parameters of the Keplerian model were leave free, including the period. The error 
bars in the RV curve are smaller than the symbols of the RV values. 
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Fig. 8.— Radial velocity curve for HR 3725, a single-line spectroscopic binary with a circular 
orbit. This plot correspond to the £t with all the parameters free, including the period. 
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Fig. 9.— Radial velocity curve for the exoplanet system HD 37605 in which we fitted all the 
parameters including the orbital period. The plot of the £t with 7 = 0 is identical to the 


eye. 















45 



0.0 0.2 0.4 0.6 0.8 1.0 

Phase 


Fig. 10.— Radial velocity curve for the exoplanet system Kepler 78b. This plot corresponds 
to the solution with circular orbit and with P and Tp fixed from the Kepler photometry. 
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Fig. 11.— Speed test results for SNR=10 synthetic datasets with different number of points. 
Each synthetic model was htted 1000 times to build 20-bin histograms. The abrupt drop for 
the lowest bin in each histogram is due to the minimum number of reannealings imposed to 
the code for termination and marks the minimum execution time for each dataset. 
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